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Abstract 

We develop a numerical method for solving a free boundary problem which de- 
scribes shape relaxation, by surface tension, of a long and thin bubble of an inviscid 
fluid trapped inside a viscous fluid in a Hele-Shaw cell. The method of solution of 
the exterior Dirichlet problem employs a classical boundary integral formulation. 
Our version of the numerical method is especially advantageous for following the 
dynamics of a very long and thin bubble, for which an asymptotic scaling theory 
has been recently developed. Because of the very large aspect ratio of the bubble, a 
direct implementation of the boundary integral algorithm would be impractical. We 
modify the algorithm by introducing a new approximation of the integrals which 
appear in the Predholm integral equation and in the integral expression for the 
normal derivative of the pressure at the bubble interface. The new approximation 
allows one to considerably reduce the number of nodes at the almost flat part of the 
bubble interface, while keeping a good accuracy. An additional benefit from the new 
approximation is in that it eliminates numerical divergence of the integral for the 
tangential derivative of the harmonic conjugate. The interface's position is advanced 
in time by using explicit node tracking, whereas the larger node spacing enables one 
to use larger time steps. The algorithm is tested on two model problems, for which 
approximate analytical solutions are available. 

Key words: Laplace's equation; Dirichlet problem; Fredholm integral equation of 
the second kind; free boundary problem, Hele Shaw flow, surface tension 



1 Introduction 



Let a bubble of low- viscosity fluid (say, air) get trapped inside a high- viscosity 
fluid (say, oil) in a quasi-two-dimensional Hele-Shaw cell: two parallel plates 
with a narrow gap between them. What will happen to the shape of the bubble, 
if the (horizontal) plates of the Hele-Shaw cell are perfectly smooth, and the 



Preprint submitted to Elsevier Science 



2 February 2008 



two fluids are immiscible? The answer depends on the initial bubble shape. A 
perfectly circular bubble (or an infinite straight strip) will not change, while a 
bubble of any other shape will undergo surface-tension-driven relaxation until 
it either becomes a perfect circle, or breaks into two or more bubbles, which 
then become perfect circles. The bubble shape relaxation process is non-local, 
as it is mediated by a flow in the external viscous fluid. The two-dimensional 
surface-tension-driven flow can be called an unforced Hele-Shaw (UHS) flow. 
This is in contrast to forced Hele-Shaw flows that have been in the focus 
of hydrodynamics and nonlinear and computational physics for the last two 
decades [1,2,3,4,5]. In rescaled variables, the UHS flow is described by the 
solution of the following free boundary problem, see e.g. Refs. [6,7]: 

V 2 p{q) = for q E E, (1) 

p(q) =K for q E 7, (2) 

v n (q) = -Vnp(g) for q E 7, (3) 

where E is an unbounded region of the plane, external to the bubble interface 
7, v n is the normal velocity of the interface, the index N denotes the component 
of vectors normal to the interface, and K is the local curvature of the interface. 
The pressure p is bounded at infinity. The free boundary problem (l)-(3) splits 
into two sub-problems: 

(1) Solving the exterior Dirichlet problem (1) and (2) and calculating v n (q) 
from Eq. (3). 

(2) Advancing the interface 7 in time with the known v n (q). 

The free boundary problem (l)-(3) represents an important example of area- 
preserving curve-shortening motion [8], but it is not integrable. Moreover, 
the only analytical solution to this problem, available until recently, was the 
approximate solution following from a linear stability analysis of a slightly 
deformed circular or flat interface [9]. Recently an asymptotic scaling theory 
has been developed for a non-trivial case when the inviscid fluid occupies, 
at t — 0, a half-infinite (or, physically, very long) strip [7]. It turned out 
that this somewhat unusual initial condition provides a useful characterization 
of the UHS flow, as the evolving strip, which develops a dumbbell shape, 
exhibits approximate self-similarity with non-trivial dynamic exponents [7]. 
Predictions of the scaling analysis have been verified numerically in Ref. [7] 
by using a boundary integral algorithm, tailored to the very large aspect ratio 
of the bubble. The present paper describes this algorithm in detail. 

A multitude of numerical methods have been suggested in the recent years for 
simulating different variants of Hele-Shaw flows. Boundary integral methods, 
which deal directly with the interface between the two fluids, are advantageous 
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compared to methods of finite elements and finite differences. Methods based 
on conformal mapping techniques have long been used in this class of problems 
(see, e.g. Refs. [2,11,12]). However, they apply most naturally to the case of 
zero surface tension and are less convenient when surface tension is non-zero 
[10]. Still another numerical strategy is phase field methods. Folch et al. [13,14] 
developed a phase field method for an arbitrary ratio of the viscosities of the 
two fluids. Unfortunately, their method becomes inefficient when the viscosity 
contrast is too high [14]. To remind the reader, the viscosity contrast is infinite 
in the case under consideration in the present work. Glasner [15] developed 
a phase field method for a description of a bubble of a high-viscosity fluid 
trapped in a low-viscosity fluid. We are unaware of any phase-field approach 
which would deal with the opposite case, which is under investigation in the 
present work: a low-viscosity bubble in a high-viscosity fluid. 

The present work suggests a numerical algorithm for solving the free bound- 
ary problem (l)-(3) in the special case of a very long bubble. It is well-known 
(but still remarkable) that the exterior Dirichlet problem (1) and (2) can be 
formulated in terms of a Fredholm integral equation of the second kind for 
an effective density of the dipole moment [16]. A naive formulation, however, 
would lead to non-existence of solution by the Fredholm alternative [17]. To 
overcome this difficulty, Greenbaum et al. [17] implemented in their algorithm 
a modification of the Fredholm equation due to Mikhlin [18]. The modified 
Fredholm equation has a unique solution for any smooth 7 and integrable /C 
[18]. Greenbaum et al. developed an efficient numerical algorithm (which is 
also valid for multiple bubbles) by discretization. However, the geometry of 
a very long and thin bubble, that we are mostly interested in, defines widely 
different length scales in the problem. Rapid variations of the dipole moment 
density at the highly curved ends of the bubble naturally necessitate a small 
spacing between the interface nodes. It is less natural, however, that, in a 
straightforward approach, one must keep the node spacing much smaller than 
the bubble thickness over the whole bubble interface. Indeed, as we show below, 
the typical length scale of the variation of the kernel of the integral equation is 
comparable to the bubble thickness which, during the most interesting part of 
the long bubble dynamics, remains almost unchanged. Apart from being com- 
putationally inefficient, the straightforward approach would cause a problem 
for explicit tracking of nodes, as the stability criterion, intrinsic in the explicit 
method, demands a time step less then a constant 0(1) multiplied by the 
node spacing cubed [19]. In this work we turned this obstacle into advantage, 
by employing the fact that the length scale of variation of the solution, over 
the most of the bubble interface, is much greater than the bubble thickness. 
We constructed a new approximation of the integral entering the Fredholm 
equation, by representing the sought dipole moment density as a piecewise 
constant function, and the bubble interface shape function as a piecewise lin- 
ear function. As a result, the integral is approximated by a sum, each term of 
which is equal to a local value of the dipole moment density multiplied by an 
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integral of the kernel between two neighboring nodes. Fortunately, the latter 
integral can be calculated analytically. The new approximation allowed us to 
considerably increase the node spacing over the most of the bubble interface, 
while keeping a good accuracy. 

Having found an approximated solution p in the form of a double layer poten- 
tial, one needs to compute the normal derivative of the solution V n p{q G 7). In 
a straightforward realization of the boundary integral formulation this would 
result in a hypersingular integral, see Ref. [17]. To overcome this difficulty, one 
resorts to theory of analytic functions and computes the harmonic conjugate 
V(q). By virtue of the Cauchy-Riemann equations, the tangential derivative 
of V(q) is equal to the desired normal derivative of p. The harmonic conjugate 
V(q) has the form of a principal value integral, over the interface, of the dipole 
moment density multiplied by a kernel, which is a function of coordinates of 
two points, q and g, belonging to the interface. This kernel diverges when the 
integration variable g coincides with q. Here we again employ the large scale 
difference at the flat part of the bubble and use the same approximation as 
in the Fredholm equation. As an additional benefit, the numerical divergence 
of the integrand of the harmonic conjugate V is avoided. As a result, we do 
not need to use even nodes to compute V at odd nodes and vice-versa, as 
suggested in Ref. [17]. 

Here is a layout of the rest of the paper. Section 2 deals with the numeri- 
cal solution of the exterior Dirichlet problem, and with the computation of 
the normal derivative of the solution at the interface. We briefly review the 
boundary integral method for an exterior Dirichlet problem and motivate the 
need for its modification for very long bubbles. Then we formulate our discrete 
approximation. In Section 3 we briefly describe a simple explicit integration 
which we used to track the bubble interface. Section 4 presents the results of 
code testing, while Section 5 presents the Conclusions. 



2 Exterior Dirichlet Problem 

2. 1 Boundary integral formulation 

Following Mikhlin [18], we seek the solution p(q) of the problem(l) and (2) 
for a simply connected bubble in a double layer potential representation: 



where /i(g) is an unknown dipole density at the point g of the interface, and 
dg is the element of arclength including the point g. The kernel K(q, g) = 




for q G E, 



(4) 
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q 



Fig. 1. Geometry of the kernel K(q,g). f and n are the tangential and outward 
normal directions, respectively. 

cosa/\f(q, g)\ follows from classical potential theory [16]. Here a is the angle 
between the outward normal to the interface at the point g and the vector 
f(q,g), see Fig. 1. 

The boundary condition (2) can be rewritten as an integral equation for /i(q), 
see, e.g., [16]: 



That is, to compute p(q) in Eq. (4), one needs to solve the integral equation 
(5). Mikhlin [18] showed that Eq. (5) has a unique solution for any smooth 
7 and integrable /C, while p(q) from Eq. (4) is a harmonic function in the 
exterior, satisfying the boundary condition Eq. (2). This representation was 
employed by Greenbaum et al. [17] for numerical analysis. 

For the purposes of the free boundary problem (l)-(3) one only needs the value 
of V n p(q), p £ 7. A straightforward calculation of V n from the double layer 
potential would yield a hypersingular integral, see below. One circumvents 
this difficulty by resorting to theory of analytic functions, see Ref. [17] and 
references therein. Suppose fj,(q) is known and introduce the quantity 



which differs from p(q) only by a constant, as §^^{g)dg = const. Obviously, 
V n p = V n p. It is known [17] that p(q) is the real part of the Cauchy integral 



-<f[l + K(q,g)]^g)d 9 = 21C(q) 

IT J-y 



(5) 
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where we have identified the points q and g on the plane with respective 



complex numbers z and (. Then p and its harmonic conjugate V satisfy the 
Cauchy-Riemann equations, so that 

Pn = V T , (6) 

where the indices n and r stand for the normal and tangential derivatives, 
respectively The kernel K(q, g) can be written as follows: 

K (q g) dg = ~ Vq) dXg + ^ 9 ~ Xq) dVa 



r 2 (q,g) 

where x and y are the Cartesian coordinates of the respective points. After a 
simple algebra we obtain 

t 7-/ \ 1 I ^{x g ,y g )(x g -x q ) u(x g ,y g )(y g -y q ) 

V{q) = ^ t r^g) ^ + r^9) ■ (7) 



2.2 Discrete approximation 



Let us parameterize the closed interface 7 of the bubble: x = x(o~), y = y(cr), 
< a < M, x(0) = x(M), y(0) = y(M), where x and y are the Cartesian 
coordinates of a point belonging to the interface. In the parametric form Eq. 
(5) becomes 

-!*(*) + - / M /i(0 [1 + «((7,0] y/± 2 +V 2 d£ = 2K{g) , (8) 

7T JO 



where 

, n _ y[x{i)-x{a)}-x[y{i)-y{a)] 
1 >*> {[x(a) - x(0\ 2 + [y(a) - yiOnV^+F ' 1 ' 

while x = dx/d^ and y = dy/dt;. The harmonic conjugate takes the form 

v w = -s7. [wo - x(°)} 2 + [wo - »wi 2 (10) 

Note that the kernel k is continuous as £ — > a. On the contrary, the integrand 
in the last expression diverges as £ — > <r, and the integral exists only as a 
principal value. In the main case of our interest the bubble length is much 
greater than its thickness A. In the almost flat parts of the interface y ~ 0. 
Now, y(o~) — y(C,) ~ A when the points o and £ belong to the different (upper 
and lower) parts of the interface, while y(o~) — ~ when they belong to 
the same part of the interface. Then, using the relation x = dx/d£, we can 
estimate the kernel k as 

-A 

[x(a) — x(^)\ 2 + A 2 
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j-1 j-1/2 j j+1/2 j+1 

Fig. 2. The discrete approximation scheme. Here £j = j. 

when a and £ belongs to the different parts of the interface, while k ~ when 
they belong to the same part. Equation (11) shows that the typical scale of 
variation of the kernel (9) over the almost flat part of the interface is of order 
of the bubble thickness A. A similar estimate applies to the fraction entering 
the integrand of Eq. (10). A straightforward discretization would then require 
a node spacing much less than A. Instead, we rewrite Eq. (8) as 

-M^) + -E/ i*m + <*,t)y& + v 2 = (12) 

where £o = 0, £ m = M, and > j = 0, 1, 2, . . . , m — 1. Introduce a 
piecewise linear approximation for and (see Fig. 2): 



(13) 



where f j < f < j = 0, 1, 2, . . . , m - 1, A;J = (x j+1 - Xj)/(^ j+1 - 

&)> b j = (Zj+iXj ~ ZjX j+1 )/(Z j+1 k) = (y j+1 - %■)/(&+! - 0). 6 i = 

(tj+iVj ~ &%-+i)/(&+i - &)» Xj = xfe), yj = yfe), and a piecewise constant 
approximation for \i: 



MO = A*j = const > < £ < £ 



(14) 



Note that < £ < ^+i) = fcj, < £ < = k). The kernel (9) is 

therefore approximated as 



where 



fc?[A££ + 6* - x(<r)] - /, ; /, /< + - 



Qj{a)/Sj 



3» = 2{k^-x(a)]+k][by y (a)]}, and C» = [6J-x(a)] 2 +[6j-y(a)] 2 . 
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The integrals in (12) can be calculated analytically: 



s 1+ 



S]e + B 3 {a)i + C 3 {a) 



1 



QM) 



2S% +1 + B J ja) 2S]j J ■ n,{a) 
arctan — —— - , x arctan — - 



2Q j (a) ■>(.), in) 

It is convenient to define £j = j, then — £j = 1. In our discretization 
scheme x(a) = x i+l/2 = kf(i + 1/2) + b x { and y(a) = y i+1 / 2 = k\(i + 1/2) + b\. 
Let us denote Qj(cr) = Qij, Bj(a) = By, and Cj(a) = Cy. The integrals in 
Eq. (12) are 



/ + K((T,0}\/i 2 +y 2 dZ = TTfij A jj , 

where 



i Jo i ^ 



Aij — — < Sj + arctan 



Kl 



vr 1 3 Q,j 4Qij + {2S]j + B, J )2Sj [ \ ■ j) + B 



uJ 



We have arrived at a set of linear algebraic equations with respect to jij, which 
is our approximation of the integral equation (5): 

m— 1 

^2 ( - 5ij)fXj = 2/Q, « = 0,1,2, 3,...,m- 1. (15) 

j=0 

We approximate the interface curvature /C(cr) by finite differences: 

K(a = . + 1/2) = JC, = [fe+1/2)2 + (i(+1/2 ) T /2 ■ 

where x i+1/2 = kf, y i+1/2 = fcf, x i+ i/ 2 = x i+2 - x m - ^ + £j-i, and j/ i+1/2 = 
Vi+2 — Vi+i — V% + Importantly, our approximation scheme yields the 
principal value of the integral (10) automatically. Furthermore, we can directly 
compute the coefficients A^-, using the same expression for i ^ j and % — j, 
where the kernel (9) has a removable discontinuity. The method suggested in 
[17] prescribes instead to use an analytic evaluation of the kernel at the point 
of removable discontinuity. 

We solved the algebraic equations (15) by an iterative refinement method after 
a LU factorization of the matrix. As the maximum number of equations in the 
examples that we considered (see below) did not exceed 1100, there was no 
need to use more sophisticated methods. 
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2.3 Grid 



Most of our results were obtained with the version of the code which assumed a 
four-fold symmetry of the bubble. This allowed us to work with a one quarter 
of the interface and reduce the number of nodes by 4. In the beginning of 
the bubble relaxation, the solution varies rapidly in the region of the lobes, 
and very slowly in the flat region of the bubble. Therefore one should employ 
here a non-uniform grid. At later times, when the aspect ratio of the bubble 
becomes comparable to unity, the code switches to a uniform grid. For the 
non-uniform grid we used an exponential spacing. Here the node spacing grows 
exponentially from the lobe's end to the middle of the flat part of the bubble. 
To generate the node distribution we use the following procedure. Let the 
quarter of the interface perimeter be IT, the specified number of nodes be m, 
and the specified smallest spacing in the lobe region be h . 

If II > ho(m — 1), the exponential grid is used. Here we introduce the quantity 
rj which satisfies the condition 

n = ho + vho + V 2 ho + ■■■+ V m ~ 2 ho = ° l - ' 1 , (16) 

I-77 

solve Eq. (16) numerically for 77, use a discrete arclength parametrization: 
£i = 0, £2 = ho, ■ ■ ■ , £k — V k ~ 2 ho, ■ ■ ■ , Cm — V m ~ 2 ho, and calculate the arrays 
x = x(£i) and y = y(£i), where i — 1, 2, . . . , m. 

In the process of the interface evolution n decreases with time, so one can 
reduce the number of the grid nodes. Furthermore, as the nodes in our code 
move like lagrangian particles (see Section 3), the node spacing in the lobe 
region decreases with time even faster. If left unattended, this would cause 
instability of the node tracking (see Section 3), as the maximum allowed time 
step is proportional to the node spacing cubed [19]. Therefore, when the min- 
imum node spacing decreases below £ 2 — O.8/10, we redistribute the nodes: we 
look for the new value of rj, corresponding to the updated value of n, cal- 
culate the new array of £, and determine the new arrays x and y by linear 
interpolation. 

When the perimeter goes down so that n < ho(m — 1), we switch to a uniform 
grid. Here we calculate a new m: m = [11/ ho] + 1, where [a] is an integer 
number such that < a — [a] < 1, and fine tune ho so that ho = H/(m — 1). 

Finally, the choice of ho is dictated by a compromise between the desired 
accuracy and the value of m which determines the size of the matrix Ay . 
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2.4 Calculation of the normal velocity 



After the set of linear equations (15) is solved, and the quantities are found, 
we compute the harmonic conjugate V. The same approximation, applied to 
Eq. (10), yields: 

j m— 1 
Vi = — 7TZ Z^i ^j^ij , 

where 

n +1 k*(k*z + bj - x t+l/2 ) + g(ge + g - yi+i/2) 

lJ I m + ^ - ^+i/2) 2 + + - 4 

1 «+i 2S 2 J + B l3 l S-?0' + l) 2 + ^-0' + l) + C ( 



2 y, 



d£=-ln 



where the quantities Sj, and C^- were defined earlier. Again, the integral 
is calculated analytically. The resulting formula for Vi is the following: 

K = -- g w In s?J . J + B(fJ . + c , (j , ("I 

where = ^(cr = i + 1/2), i = 0, 1, 2, 3, . . . , m - 1, see Fig. 2. Note that for 
£ = i + 1/2 the denominator of the integrand in Fij vanishes, and the integrand 
diverges. To overcome this problem, Ref. [17] suggested to divide the mesh 
into odd and even nodes and compute V at the odd points by summing over 
the even nodes, and vice-versa. Our analytical integration yields the correct 
principal value of the integral, so there is no need to use the recipe of Ref. 
[17]. " 

To compute the normal velocity of the interface we use the Cauchy-Riemann 
equation (6) and approximate the derivative of V with respect to the arclength: 

v n {a = i) = z , 



where 

yj (x i+1 - Xi) 2 + (y i+ i - ViY + \J (xi - Xi_i) 2 + (yi - 2/i-i) 2 



1 

Si — - 
2 



3 Interface Tracking 

To track the interface, we use an explicit first-order integration: 

Xi(t + At) = Xi(t) + At v n (a = i, t) cosnj, 
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yi(t + At) = yi(t) + Atv n (a = i,t) sinn, 



where 



cosn, = 




smrij = 




Xi = (xj+i/2 + Xi_i/ 2 )/2 = (fcf + kf_ 1 )/2, and = (fcf + fcf_i)/2. We have 



assumed the counter-clockwise direction of the interface parametrization, see 



It is important to prescribe the time step At properly. We employ an ad-hoc 
criterion which demands that the node displacement at each grid point be con- 
siderably less then the curvature radius Ri at that point: min | (At v n (i))/Ri\ < 
e, < i < m — 1. That is, we consider the curvature radius Ri as a natural 
local length scale of the problem. A more convenient form of this criterion is 



where e is an input parameter which has to be sufficiently small to satisfy 
the requirements of stability and desired accuracy. In the exact formulation 
(l)-(3) the bubble area must be constant in the process of relaxation. The 
area conservation can be conveniently used for accuracy control of the code. 



4 Numerical Results 

We present here some simulation results produced with our code for two dif- 
ferent sets of initial conditions. One of them describes the decay of a small 
sinusoidal perturbation of a perfectly circular bubble of inviscid fluid. An ap- 
proximate analytical solution to this problem is given by the linear stability 
analysis [9], and we used this solution to test the code. 

The second initial condition describes a very long and thin strip of inviscid 
fluid. In the process of its shrinking the bubble develops a dumbbell shape, 
while the characteristic dimensions of the dumbbell exhibit asymptotic scaling 
laws found in Ref. [7]. 

4-1 Relaxation of a slightly perturbed circle 

Let the initial shape of the interface be a circle with a small sinusoidal per- 
turbation: 



where p and ip are the polar radius and angle, respectively, Rq is the radius 
of the unperturbed interface, while 5(0) and n are the initial amplitude and 



Fig. 1. 



At = emin{\Ri/v n (i)\} , 




p{<p, 0) = Ro[l + 5(0) sin(ny9)], 
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Fig. 3. Shown in the logarithmic scale is the perturbation amplitude 5 as a func- 
tion of time. The squares are the simulation results, the solid line is the analytical 
prediction. 

azimuthal number of the perturbation. The analytical solution provided by 
the linear theory [20] is 



p(ip, t) = Rq[1 + 6(t) sin(ny?)] , 



where the amplitude of the perturbation is 

n(n 2 — 1) 



S(t) = 5(0) exp 



A typical numerical result is presented in Fig. 3. The parameters are Rq = 100, 
5(0) = 0.01 and n — 4. In the case of n = 4 the interface has a four-fold 
symmetry which allows a direct application of our code. In this simulation the 
quarter of the interface was described by 100 nodes. The initial spacing was 
uniform. The code did not have to use the mesh interpolation in this example. 
The parameter regulating the time step was e — 5 • 10~ 5 . As one can see, a 
very good agreement with the analytical result is obtained. 



4-2 Relaxation of a long and thin bubble 



In the second setting the initial interface shape is a very long rectangular strip. 
In the example we report here the initial strip thickness was 1, and the initial 
length 2000. Here we could compare the numerical results with the predictions 
of a recent asymptotic scaling analysis [7]. The interface shapes at different 
times are presented in Fig. 4. It can be seen that the shrinking strip acquires 
the shape of a dumbbell (or petal). At much later times it approaches circular 
shape. By the end of the simulation (at t = 48000) the relative deviation of the 
observed shape from the perfect circle, [p max ((p) — p m i n (v)}/ Pmin(v) ~ 0.013. 
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200 400 600 800 1000 455 460 465 470 475 20 40 60 



x xx 

Fig. 4. Figure a shows a snapshot of one half of the simulated system at t = 0, 3670, 
7020, and 24840. Notice the large difference between the horizontal and vertical 
scales. Figure b shows the lobe of the dumbbell to scale at t = 7020, while Figure c 
shows the computed bubble shape at late times: t = 30900, 34200 and 48000. 



100: 




10 I 1 I ■ — 4 

100 1000 100 1000 
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Fig. 5. Figure a shows, in a log-log scale, the retreat distance L versus time and its 
power- law fit 2.75 1 60 . Figure b shows, in a log-log scale, the maximum dumbbell 
height, h max (the empty circles), and the position of the maximum, x™ ax (the 
filled circles), versus time, as well as their power-law fits 0.66 i - 21 and 0.94 1 - 20 , 
respectively. 

The asymptotic scaling analysis [7] deals with the intermediate stage of the 
relaxation. Introduce the retreat distance L(t) = 1000 — xu p (t) : where x t i P (t) is 
the maximum abscissa of all points belonging to the interface. One prediction 
of Ref. [7] is that, at intermediate times, L(t) oc t 3 ^ 5 . Figure 5a shows a 
very good agreement of this prediction with the simulation result. Additional 
predictions of asymptotic scaling analysis deal with the time dependence of the 
maximum dumbbell elevation h max (t), and of the abscissa of the corresponding 
point of the interface x max (t). Let us introduce a new variable: x±(x, t) = 
Xu P {t) — x, the distance along the x-axis between the tip of the dumbbell 
and a point x. In particular, x™ ax {t) = x t i P (t) — x max (t). A comparison of the 
simulation results with the predicted intermediate-time scaling laws h max (t) oc 
x™ ax (t) oc t 1 / 5 is shown in Figure 5b, and again a very good agreement is 
observed. 

To verify the self-similarity of the dumbbell shape in the lobe region, predicted 
in Ref. [7], we introduce a new function h(x 1 ,t) so that h[xi(x,t),t] = y(x,t). 
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Fig. 6. Self-similarity of the lobe. Shown is the shape function h{x\,t), rescaled to 
the maximum dumbbell elevation, versus the coordinate x%, rescaled to the abscissa 
of the maximum, at times 160.3 (the filled circles), 1000 (the squares), and 3010 
(the empty circles). 

Figure 6 shows the spatial profiles of h rescaled to the values of h max versus 
Xi/x™ ax at three different times. The observed collapse in the lobe region 
confirms the expected self-similarity. 

The initial number of nodes in this simulation was 1100, and the smallest 
spacing in the lobe region was 0.4. With the grid interpolation employed, the 
time-step parameter e = 0.005 proved sufficiently small to guarantee stability 
and good accuracy. As the curvature of the interface goes down during the 
evolution, the required time step increases significantly. It was 1.7 x 10~ 3 at 
t = 0, 0.22 at t = 3670 and increased up to about 10 by the end of the 
simulation, at t — 48000. We used the small observed area loss of the bubble 
for accuracy control. The observed area loss was less then 0.5% for t < 10000. 
By the end of the simulation, at t — 48000, the area loss reached only 2.8%. 
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5 Conclusion 



We have developed and tested a new numerical version of the boundary in- 
tegral method for an exterior Dirichlet problem, which is especially suitable 
for long and thin domains. The method allows one to significantly reduce the 
number of the interfacial nodes. The new method was successfully tested in a 
numerical investigation of the shape relaxation, by surface tension, of a long 
and thin bubble, filled with an inviscid fluid and immersed in a viscous fluid 
in a Hele-Shaw cell. Here we confirmed the recent theoretical predictions on 
the self-similarity and dynamic scaling behavior during an intermediate stage 
of the bubble dynamics. 
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